1. Field of the Invention
Embodiments of the present invention relate generally to methods of simulating deformable objects. More particularly, embodiments of the invention relate to methods of simulating deformable objects using a geometrically motivated underlying model.
2. Description of Related Art
Realistic computer simulations of everyday objects such as clothing, plastics, elastic materials, and articulated joints can be extremely difficult to generate due to the complex ways in which these objects tend to deform in real life. Simulations of such “deformable objects” must generally take into account both complicated geometric features and material properties of the objects.
One common technique used to simulate deformable objects first creates a virtual model of an object (e.g., a mesh or point cloud) and then applies simulated physical forces such as tension, friction, gravity, pressure, etc., to the discrete points of the virtual model. Such virtual models have been used to represent a wide variety of materials under different conditions. For example, researchers have developed virtual models for clothing, plastic, rubber, and so on. In addition, researchers have also developed virtual models to simulate complex, unique behaviors of these objects, such as fracturing and melting.
Some of the more common approaches to simulating deformable objects involve finite difference methods, mass-spring systems, boundary element methods, finite element methods, and implicit surfaces and mesh-free particle systems.
Despite the number of methods that have been developed for simulating deformable objects, such methods are rarely incorporated into applications like computer games. For example, a few contemporary games incorporate cloth models with simple geometries; however, in general, most applications running on Personal Computers (PCs) or game consoles tend to model only rigid bodies or objects. Indeed, multiple rigid objects are sometimes combined to imitate the movement of a deformable object, but true deformable object simulations are rare.
One reason why deformable object simulations are rarely used in an application such as computer games is that stable simulations are generally too inefficient to satisfy the real-time demands of the application. For instance, conventional simulation methods often rely on implicit numerical integration to update the locations of discrete points in a virtual model. Implicit numerical integration techniques can provide a stable simulation, but they are too computationally expensive to process large and complicated physical models in realtime. Other physical modeling approaches rely on explicit numerical integration, which is more efficient, but it cannot guarantee a stable simulation.
The term “stability” here refers to a process's tendency to respond in a reasonable way to minor deviations in its inputs. For instance, implicit numerical integration techniques produce accurate simulations as various input parameters are varied, such as the mass of simulated points, the timestep of the integration, and so on. In contrast, explicit numerical integration techniques may produce simulations where the overall energy of a system erroneously increases by simply varying the timestep of the integration, the mass of simulated points, or the stiffness of a simulated object. As a result, explicit numerical integration techniques can produce highly unrealistic simulation results.
Various approaches have been developed to improve the efficiency of stable simulation techniques. For example, robust integration schemes with large timesteps and multi-resolution models have been developed. In addition, modal analysis approaches, which trade accuracy for efficiency, have also been developed. Furthermore, methods incorporating pre-computed state space dynamics and precomputed impulse response functions have also been used to improve the efficiency of simulations. Finally, dynamic models derived from global geometric deformations of solid primitives such as spheres, cylinders, cones, or super quadrics have also been introduced to improve the efficiency of stable simulation techniques.
FIG. 1 is a diagram illustrating one way in which instability arises in deformable object simulations using explicit numerical integration. In FIG. 1, a simple one dimensional deformable object is modeled as a mass-spring system. Like many physically motivated deformable object simulations, the mass-spring system relies on Newton's second law of motion. According to Newton's second law of motion, the acceleration of an object produced by a net force is directly proportional to the magnitude of the net force in the direction of the net force, and inversely proportional to the mass of the object. In deformable objects, at least part of the net force at any point is created by displacement of the point from an equilibrium position. The displacement of a point from its equilibrium position creates potential energy, i.e., “deformation energy”, causing the point to be pulled toward the equilibrium position.
Referring to FIG. 1A, a mass-spring system 100 comprises a spring 101 with a resting length l0, and two point masses 102 and 103 both having mass “m”. Point mass 102 is fixed at an origin and point mass 103 is located at x(t) at an initial time “t”. At initial time “t”, the amount of force on point mass 103 is defined by the spring equation f=−k(x(t)−l0), where “k” is the spring constant, or “stiffness”, of spring 101. Conceptually, the spring equation indicates that point mass 103 is pulled toward an equilibrium position where x=l0.
The location of point mass 103 is updated by a modified Euler integration scheme after a timestep “h”. According to the modified Euler integration scheme, the velocity “v” and position “x” of point mass 103 at time “t+h” are computed using the following equations (1) and (2):
                              v          ⁡                      (                          t              +              h                        )                          =                              v            ⁡                          (              t              )                                +                      h            ⁢                                          -                                  k                  ⁡                                      (                                                                  x                        ⁡                                                  (                          t                          )                                                                    -                                              l                        0                                                              )                                                              m                                                          (        1        )                                          x          ⁡                      (                          t              +              h                        )                          =                              x            ⁡                          (              t              )                                +                      hv            ⁡                          (                              t                +                h                            )                                                          (        2        )            Equation (1) uses an explicit Euler step and equation (2) uses an implicit Euler step.
Equations (1) and (2) can be represented as a system matrix “M” multiplied by a state vector [v(t),x(t)]T, i.e.,
            [                                                  v              ⁡                              (                                  t                  +                  h                                )                                                                                        x              ⁡                              (                                  t                  +                  h                                )                                                        ]        =                  [                                                            v                ⁡                                  (                  t                  )                                                                                                        x                ⁡                                  (                  t                  )                                                                    ]            +                        1          m                ⁡                  [                                                                                          kl                    0                                    ⁢                                      h                    2                                                                                                                                            kl                    0                                    ⁢                  h                                                              ]                      ,where system matrix “E” is defined by the following equation (3):
                    E        =                  [                                                    1                                                              -                                      kh                    m                                                                                                      h                                                              1                  -                                                                                    h                        2                                            ⁢                      k                                        m                                                                                ]                                    (        3        )            System matrix “E” has eigenvalues e0 and e1 represented by the following equations (4) and (5):
                              e          0                =                  1          -                                    1                              2                ⁢                                                                  ⁢                m                                      ⁢                          (                                                                    h                    2                                    ⁢                  k                                -                                                                                                    -                        4                                            ⁢                                                                                          ⁢                                              mh                        2                                            ⁢                      k                                        +                                                                  h                        4                                            ⁢                                              k                        2                                                                                                        )                                                          (        4        )                                          e          1                =                  1          -                                    1                              2                ⁢                                                                  ⁢                m                                      ⁢                          (                                                                    h                    2                                    ⁢                  k                                +                                                                                                    -                        4                                            ⁢                                                                                          ⁢                                              mh                        2                                            ⁢                      k                                        +                                                                  h                        4                                            ⁢                                              k                        2                                                                                                        )                                                          (        5        )            Since system matrix “E” represents a discrete system, the spectral radius of system matrix “E”, i.e., the maximum magnitude of eigenvalues e0 and e1, must not be larger than one (1) to ensure stability of the discrete system. The magnitude of eigenvalue e0 converges to 1 with |e0|<1 for h2k→∞. However, the magnitude of e1 is only smaller than one where timestep “h” is less than
  2  ⁢                    m        k              .  Where timestep “h” is greater than
      2    ⁢                  m        k              ,the system is unstable. Accordingly, the integration scheme involving equations (1) and (2) is only conditionally stable.
To further illustrate the instability of the discrete system represented by system matrix “E”, FIG. 1B shows the result of performing an integration step starting with v(t)=0. The integration step moves point mass 103 by a distance
      Δ    ⁢                  ⁢    x    =            -                                    h            2                    ⁢          k                m              ⁢                  (                              x            ⁡                          (              t              )                                -                      l            0                          )            .      Where timestep “h” or stiffness “k” are too large or mass “m” is too small, point mass 103 overshoots equilibrium position l0, by a distance greater than the distance between x(t) and l0. In other words, |x(t+h)−l0|>|x(t)−l0|. As a result, the potential energy of system 100 increases after timestep “h”. Since system 100 had zero kinetic energy at time “t”, the overall (i.e., kinetic plus potential) energy of the system is erroneously increased after timestep “h”.
In general, the stability problem of explicit integration schemes can be stated as follows: elastic forces are negative gradients of elastic energy. As such, elastic forces point towards equilibrium positions. Explicit schemes may inaccurately scale the elastic forces when computing the displacements of points, causing the points to overshoot equilibrium positions so much that they increase the deformation and the energy of the system instead of preserving or decreasing the deformation energy, which is required for stability.
One way to address the overshoot problem is to limit the displacement of points such that they never overshoot their respective equilibrium positions. For instance, in the one-dimensional spring example of FIG. 1, the movement of point mass 103 could be restricted from passing the equilibrium position x=l0. One problem with this approach is that for many types of physical models, equilibrium positions are not readily defined for all points. For example, it is difficult to define equilibrium positions in solid finite elements or geometrically complex meshes.